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Abstract 

It is often useful to represent the infectious dynamics of mobile agents by metapopu- 
lation models. In such a model, metapopulations form a static network, and individuals 
migrate from one metapopulation to another. It is known that heterogeneous degree 
distributions of metapopulation networks decrease the epidemic threshold above which 
epidemic spreads can occur. We investigate the combined effect of heterogeneous degree 
distributions and diffusion on epidemics in metapopulation networks. We show that for 
arbitrary heterogeneous networks, diffusion suppresses epidemics in the sense of an in- 
crease in the epidemic threshold. On the other hand, some diffusion rates are needed to 
elicit epidemic spreads on a global scale. As a result of these opposing effects of diffusion, 
epidemic spreading near the epidemic threshold is the most pronounced at an intermedi- 
ate diffusion rate. The result that diffusion can suppress epidemics contrasts with that 
for diffusive SIS dynamics and its variants when individuals are fixed at nodes on static 
networks. 
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1 Introduction 



Infectious diseases are transmitted through social contacts between individuals. The relation- 
ships between the structure of social networks and the number of infected individuals have 
been extensively studied in mathematical epidemiology and network sciences. Such studies 
have established that heterogeneity in contact rates enhances epidemics [IHl]- 

The contact networks underlying, for example, sexually transmitted diseases of humans 
and computer viruses on the Internet are considered to be static on the time scale of epi- 
demics. However, the social networks underlying humans' prevailing infectious diseases such as 
influenza are presumably dynamic even during a day. The dynamics of networks are induced 
by the migration of individuals among residences, workplaces, places for social activities, and 
so on. Other animals also migrate between habitats. Metapopulation models are useful in 
describing epidemics and ecological invasions in such a situation [Il|2l[5l[6]. A node in such a 
model represents a metapopulation or a habitat, and not an individual. A link represents a 
physical pathway connecting a pair of metapopulations. Individuals travel from one node to 
another. Simultaneously, interactions between individuals, such as infection, can occur in each 
metapopulation. 

The basic reproduction number and the condition for the occurrence of global epidemics 
have been theoretically determined for the susceptible-infected-suscepetible (SIS) dynamics 
and its variants on metapopulation networks with general connectivity profiles [2|[7]-[T0]. Real 
metapopulation networks relevant to epidemics, such as networks of urban metapopulations 
and networks of airports, often have heterogeneous degree distributions; some metapopulations 
are highly connected as compared to many others [SHUITT] . Colizza et al. put important efforts 
for understanding infectious dynamics on complex metapopulation networks . In partic- 

ular, they showed that heterogeneous degree distributions enhance epidemics in uncorrelated 
networks of metapopulations (also see [IT]). Similar results have been established for models 
of epidemics on static networks of individuals (see [2111] for reviews). However, the relationship 
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between the epidemic threshold of the infection rate (or population density) above which large 
epidemics can occur and the degree distribution differs in these two classes of models. 

For metapopulation models, effects of stochasticity in infection, recovery, and migration dy- 
namics were also analyzed in a recent paper [18] . Other authors also investigated the threshold 
for infection and epidemic dynamics in uncorrelated heterogeneous networks of metapopula- 
tions [UlEO]. 

In the present work, we investigate the effect of diffusion on epidemics in metapopulation 
networks on which individuals diffuse. Diffusion increases the possibility of epidemics for the 
SIS model [2T1|22] and its variants [23|l2^ in conventional static networks of individuals. This 
is presumably because diffusion helps infected individuals to contact new susceptibles to be 
infected. The effect of traffic-driven diffusion on the epidemic threshold was also studied recently 
[25] . We show that for epidemics in metapopulation networks, diffusion prohibits rather than 
enhances epidemic spreads in the sense of the epidemic threshold. Although this result has 
been implicitly recognized in previous literature [20], we show it by developing an analytical 
method to calculate the epidemic threshold for arbitrary networks and diffusion rates. We then 
support this result by numerical simulations and further show that an intermediate diffusion rate 
magnifies epidemics near the epidemic threshold. Qualitatively identical behavior is observed in 
numerical simulations of the susceptible-infected-recovered (SIR) dynamics on metapopulation 
networks. 

2 Model 

We assume connected and undirected networks of metapopulations; extending the following 
results to the case of weighted networks is straightforward. Each node represents a metapopu- 
lation or habitat and accommodates any number of individuals. For a network having nodes, 
the N hj N adjacency matrix is denoted by A; Aij = 1 when node i and node j are adjacent, 
and Aij = otherwise. A is symmetric and we set An = {1 < i < N). 

We consider the SIS dynamics of diffusive individuals on this network using the continuous- 
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time formulation developed in [T^[2D]. The use of the original discrete-time formulation of 
the model [13] does not essentially change the following results. The network contains N p 
individuals that independently diffuse on the network. Each individual stays at a node and 
takes either a susceptible or an infectious state. The SIS dynamics with infection rate /3 and 
recovery rate /i occur in each metapopulation. We assume mass interaction and therefore do not 
normalize the infection rate by the number of individuals in the metapopulation [13] . By doing 
so, we focus on situations in which metapopulations that host relatively many individuals would 
be a source of epidemics once an infected individual is introduced to such a metapopulation. 
The infection event at node i occurs at a rate of (3ps.iPi,i, where ps.i and are the numbers of 
susceptible and infected individuals at node i, respectively. This assumption implies all-to-all 
interaction within each metapopulation. 

At the same time, individuals perform a simple random walk. In infinitesimal time At, a 
susceptible (infected) individual at node i with degree ki moves to one of its neighboring nodes 
with equal probability ~ DsAt/ki DiAt/ki), where Ds (-Di) is the diffusion rate for the 
susceptible (infected) individual and ki is the degree of node i. 



3 Epidemic threshold for arbitrary metapopulation net- 



We derive the epidemic threshold above which the endemic state {i.e., survival of infected 
individuals) can occur. Although general solutions have been formulated [2|[71-[T0] and solutions 
are known for uncorrelated random networks [T3lfT8] - [20] . we are concerned with the effect of 
the diffusion rate on the epidemic threshold in general heterogeneous networks. 
The master equations are given by 



works 



dps,i 
dt 



(3ps,iPi,i + m,i - Dsps,i + I^s -^Ps 



(1) 




(2) 
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To derive the epidemic threshold, we consider the steady state in which pi i (1 < i < A^) is in- 
finitesimally small. In this situation, Eq. ([T]) represents the master equation for the continuous- 
time pure diffusion of susceptible individuals because Pps,iPi,i ~ and yupi j ~ 0. Given > 0, 
the steady state for the number of susceptible individuals is given by 

Ph = (3) 

where (k) = ^i/^ is the average degree and p = ^^i(ps,j + Pi,i)/N, the total population 
density. Substituting Eq. in Eq. yields 

^P^'i f^P 1 n , n ^Ji fA\ 

\ / j=i i 

Eq. (jlj) has Pi j = (1 < ? < A^) as the disease-free solution. The destabilization of the disease- 
free solution implies an endemic state, that is, pj^ > (1 < z < A^) [T1[T0]. We focus on the 
threshold value of /3 above which the endemic state is induced, which is denoted by /3c. We note 
that the threshold obtained below can also be expressed in terms of p, as is done in previous 
literature [l3l[T9l[20] , because (3 and p appear only as the multiple /3p in Eq. (j4]). 
In the strong diffusion limit = oo, Eq. (jl]) implies 

k- 

Pi,i = j^Ph (5) 

where pi = J2f=iPi,i/-^ ~ is the total density of the infected individual. By substituting 
Eq. (|5]) in Eq. (jl]) and taking the summation over i, we obtain 

It = [-IJ^ - («> 

which reproduces the threshold (3^ = p (k)'^ / p {k"^) known for uncorrelated networks [13]. When 
Di = oo, heterogeneous degree distributions decrease the epidemic threshold in arbitrary 
metapopulation networks. 

When Di = 0, Eq. @ represents decoupled dynamics, and /3c is equal to the epidemic 
threshold at an isolated node i that hosts the largest number of individuals among the A^ nodes. 
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We obtain 

A = (7) 

P"'max 

where kma.x is the degree of node i [20]. When (3 > f3c, the nodes that satisfy ki > (k) fi/{f3p) 
can eventually have infected individuals, whereas the other nodes can not. In heterogeneous 
networks, we obtain (k^) / {k) = Y^iLi kf / Y!j=i kj < YliLi kikms.^/ Y.j^=i kj = /cmax- Therefore, 
the value of (3c for Z^i = is smaller than that of /3c for Di = oo. Although this fact apparently 
indicates that the endemic state occurs more easily for Di = than for Di = oo, this state- 
ment is somewhat misleading because Di = implies that the infection does not travel across 
metapopulations. When Di = 0, the infection is confined in the metapopulations that contain 
index patients. 

Assume that Ds and Di are arbitrary. We develop a method to calculate the epidemic 
threshold for an arbitrary structure of networks and diffusion rates. The Jacobian J = {Jij) of 
the dynamics represented by Eq. (jlj) around the disease-free solution is given by 

J.,= (^-/i-A)5., + A^, (i<^,j<iV) (8) 

where 5ij is the Kronecker delta. When the real parts of all the eigenvalues of J are negative, 
the disease-free solution is stable. Otherwise, the endemic state emerges. Similar treatments 
were carried out for this model on uncorrelated random networks [20] and the SIS dynamics on 
static networks of individuals [261 - 129] . 
J is isospectral to 

J' = ( JM = diag Jdiag f v^, . . . , , (9) 



where diag((ii, . . . , djv) denotes the diagonal matrix with the zth diagonal element equal to d 
Because 

is symmetric, all the eigenvalues of J' and hence those of J are real. Let Amax denote the 
maximum eigenvalue of J. 
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J'.^=(^-f^-Di)s.,+D,^ (10) 



Fix (3 and express Amax using the Rayleigh quotient: 

Amax = maxx^ J'x, (x E R^) (11) 

\x\ = l 

where T denotes the transpose. Suppose that the maximum of Eq. flTTl) is reahzed hj x = x, 
where \x\ = 1. When /3 is replaced by /3 + A/3, we obtain x^ J'x = Amax + pA/3 X^ili ^j^f/ 
where Amax is the largest eigenvalue before (3 is replaced by /3 + A/3. Therefore, Amax niono- 
tonically increases with /3; this guarantees that the minimum value of /3 satisfying Amax = is 
equal to (3c- The condition Amax = is equivalent to 

detJ = 0. (12) 

With regard to the value of (3c, Eq. f lT^ is equivalent to 



0=Wdet 

P 



Equation (fT3|) indicates that the minimum eigenvalue of the symmetric matrix M = (Mij), 
where 

fi + Di DiAij\ (k) 



is equal to /3c. For arbitrary networks, we can reliably compute the minimum eigenvalue of 
M in 0{N^) time by combining the Cholesky decomposition, which is applicable to symmetric 
matrices, and the inverse power method. 

For uncorrelated networks, substituting Aij = kikj/{{k) N) in Eq. (|T4|) yields 

The inverse of Eq. (fTSi) is given by 

Because all the eigenvalues of M are positive, all the eigenvalues of are also positive. 
Therefore, we can rapidly calculate (3c by applying the standard power method to M~^. 



4 Epidemic threshold increases with the diffusion rate 

In this section, we show that /3c monotonically increases with Di in arbitrary heterogeneous 
networks. To this end, we decompose M as 

M={B + DjL)^, (17) 



where 

B = nx diag ( , . . . , 
and 

1 1 \ 



5 = //xdiagf-^,...,-^ ) (18) 



(y/ci ' ■ ■ ■ ' /cat J kikj ' ^^^"^ 
L is a Laplacian matrix such that its minimum eigenvalue is 0. All the other eigenvalues of 
L are positive for connected networks. The interlacing eigenvalue theorem (also called Weyl's 
inequality) implies that the minimum eigenvalue of the summation of two symmetric matrices 
is not smaller than the sum of the minimum eigenvalues of the two individual matrices [5U|IHT]. 
Because the minimum eigenvalue of B is equal to /S^ for Z^i = and that of L is equal to 0, the 
minimum eigenvalue of M, i.e., for a given Di > 0, is not smaller than for Di = 0. More 
generally, if D[> Di, we can apply the abovementioned arguments to the sum of -B + DiL and 
{D[ — Di)L, i.e., B + D[L. Therefore, /3c monotonically increases with Di. 

For regular graphs, that is, networks in which all the nodes have the same degree, Pc = p/p 
both for Di = and Z^i = oo. Because /3c monotonically increases with Di, 13^ = fi/p holds 
true for any Di > 0. The diffusion rate affects the epidemic threshold for heterogeneous 
metapopulation networks but not for homogeneous networks. 

5 Numerical results 

The results presented in the previous section imply that diffusion suppresses epidemics on het- 
erogeneous metapopulation networks in the sense of an increased epidemic threshold. However, 
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if the diffusion rate is too small, infected individuals may be localized within the metapopula- 
tions having index patients. We investigate the combined effect of the heterogeneity and the 
diffusion rate by carrying out direct numerical simulations of infectious dynamics. 

We set the recovery rate /i = 1 and population density p = 50 without loss of generality 
and vary the infection rate /3. For simplicity, we set D = Ds = Di. We start with one 
infected individual selected from the population with equal probability 1/Np. The other Np — 1 
individuals are initially susceptible. We measure the fraction of infected individuals in the 
steady state, denoted by pj. We calculate pj as the number of infected individuals after transient 
divided by and averaged over multiple realizations. The relaxation time is governed hj 1/D 
for small D (> 0) and 1/p otherwise. Therefore, we set the transient length to 500 for = 
and > 4, and 2000/L' for < D < 4. 

We first examine the SIS dynamics on scale-free networks generated by the Barabasi-Albert 
(BA) model [32] with A^ = 200 and mean degree 6 [i.e., 3 links are added per new node). To 
generate a network, A^ — 3 = 197 nodes are sequentially added to a triangle according to the 
preferential attachment. The degree distribution of the network is the power law p{k) oc 
|32] . For a generated scale- free network, in Fig. [T](a), we compare /3c obtained theoretically and 
numerically for a range of D. The theoretical value of Pc is equal to the minimum eigenvalue 
of matrix M given by Eq. (I14p . We obtain a numerical estimate of /3c as the value of /3 above 
which an infected individual survives in at least one of the 5000 realizations. The numerical 
results are in a good agreement with the theoretical results. 

Next, we examine the SIS dynamics for suprathreshold /3. For fixed /3 and D, we generate 
400 realizations of the scale-free network and run the SIS dynamics 5 times. The mean fraction 
of infected individuals pj is calculated on the basis of 2000 runs. The dependence of pj, 
normalized by p, on /3 and D are shown in Fig. [It^b). In agreement with Fig. [I](a), the epidemic 
threshold increases with the diffusion rate. 

However, when D = 0, pi is small for any /3. This is because infected individuals do not 
migrate to other metapopulations, and they are localized within the single metapopulation 
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to which the index patient belongs even above the epidemic threshold. When D = 0.1, the 
infection reaches multiple metapopulations above the epidemic threshold, but only to a limited 
extent because of a small diffusion rate. When D = 1, above the epidemic threshold, the 
infection seems to spread to more metapopulations and individuals than when D = 0.1. A 
visible fraction of infected individuals (pj/p ~ 0.03) survives between /3 ^ 0.4 and (3 ^ 0.6, 
for which larger diffusion rates do not allow the endemic state. When D = 20, pj is large 
for large values of f3, presumably because diffusion helps the infected individuals to reach 
metapopulations with small degrees. We stress that the epidemic threshold for = 20 is 
nonetheless larger than that for D = 0.1 and D = 1, which is consistent with Fig. [^a). 

To quantify the localization of infected individuals, in particular when an appropriately 
small D [i.e., D = 1) yields relatively many infected individuals {i.e., 0.4 < /3 < 0.6), we 
plot the inverse participation ratio defined by Y2 = Yld=i{Pi,i/ PiY (see [33] for an exemplary 
use of I2 for epidemic dynamics on networks). When Y2 is of the order of unity, the infected 
individuals are localized in a small number of metapopulations. When Y2 = 0{1/N), the 
infected individuals have spread to 0{N) metapopulations. In Fig. HJ^c), Y2 are plotted for the 
values of D used in Fig. Hl^b). For D = 0, Y2 = 1 irrespective of /3, as expected. For D = 0.1, 
Y2 is relatively large for /3 ^ 0.5. For D = 1, Y2 in this range of /3 becomes smaller to approach 
the values for D = 20 (note the logarithmic scale in Fig. [II^c)). 

In infectious diseases with which metapopulation models are concerned, such as influenza, 
the SIR model seems to have a wider applicability than the SIS model [T21[I11I31]- In the SIR 
model, the individuals that recover from the infected state do not return to the susceptible 
state but transit to the recovered state. When the basic reproduction number is assumed to 
be the same for each metapopulation, analytical results for the SIR model on heterogeneous 
metapopulation networks have been established [T5|ll6pi8] . However, the analysis seems difficult 
when the basic reproduction number depends on the metapopulation. In our model, highly 
populated metapopulations, which presumably make contact with many other metapopulations, 
have large basic reproduction numbers. 
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We carry out numerical simulations starting from an arbitrarily chosen infected individual 
until there is no infected individual. Then, we measure the final fraction of the recovered 
individuals averaged over 2000 realizations. The results shown in Fig. are qualitatively 
the same as those for the SIS model shown in Fig. [T](b). To confirm that an appositely small D 
facilitates spreads of infection to many metapopulations, we measure Y2. For the SIR dynamics, 
we need to modify the definition of I2 because runs without any secondary infections lead to 
Y2 = 1. Runs with a minor fraction of infected individuals also yield a spuriously large Y2. 
Therefore, in the calculation of Y2, we exclude the runs in which less than 1% of the individuals 
are eventually infected. The dependence of Y2 on /3 and D shown in Fig. Et^b) are qualitatively 
the same as those for the SIS model shown in Fig. HJ^c). 

As an example of real heterogeneous networks, we simulate SIS and SIR dynamics on the 
network of US airports yjj. This network has = 332 nodes, 2126 links, and a long-tailed 
degree distribution. We ignore the link weight. A airport is considered to be a metapopulation 
that hosts traveling individuals. Figure [3I^a) indicates that the values of Pc obtained from direct 
numerical simulations are in a good agreement with the theoretical values. The suprathreshold 
behavior of the SIS and SIR dynamics shown in Figs. ElJ^b) andEI^c), respectively, is qualitatively 
the same as that on the BA model (Figs. [Hb) and[2]^a)). 

The combined effect of the heterogeneity and the diffusion rate is absent for homogeneous 
networks. To demonstrate this, we carry out numerical simulations on the regular random 
graph with A^ = 200 and degree 6; all the nodes have degree 6 and are wired randomly under 
the condition that the generated network is connected. The theory obtained in Sec. H] indicates 
that the epidemic threshold Pc is independent of D for the regular random graph. As shown in 
Fig. m^a), the numerically obtained /3c is largely independent of D and close to the theoretical 
value, albeit there is some disagreement for small D. For a range of /3, the normalized pj and 

for the SIS and SIR models are shown in Figs. ID^b) andlD^c), respectively. Because /3c little 
depends on D, a large diffusion rate results in a large fraction of infected individuals for any /3. 
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6 Discussion and Conclusions 



We have shown that diffusion increases the epidemic threshold of the SIS dynamics in arbitrary 
heterogeneous networks. Nevertheless, a certain amount of diffusion is necessary to enable the 
transmission of epidemics across metapopulations. Numerical simulations of the SIR dynamics 
also yielded similar results. Similar conclusions are expected for other population dynamics 
used in ecology [B] . Indeed, dispersal decreases the population growth rate in ecological invasion 
models in which each metapopulation carries heterogeneous birth and death rates [55] . 

Our numerical results indicate that excessive diffusion suppresses epidemics in metapopu- 
lation networks near the epidemic threshold. This effect of diffusion contrasts with that on the 
diffusive SIS model and variants in static networks of individuals [2T1 - [2l] . We note that, in 
the metapopulation model, an increased diffusion rate can also suppress the spreading speed 
of epidemics. In heterogeneous static networks of individuals, the spreading speed in an early 
stage of transmission is controlled by the largest eigenvalue of the matrix that represents trans- 
mission of infection |[33j. The counterpart in the metapopulation model is the matrix given 
by Eq. ([8]). The largest eigenvalue of the matrix given by Eq. (|8]) decreases with the diffusion 
rate; this can be shown by adapting the proof of the monotonic dependence of the epidemic 
threshold on the diffusion rate (Sec. H]). 

The nonmonotonic dependence of the fraction of infected individuals on the diffusion rate 
near the epidemic threshold is derived from the fact that epidemics confined in a single metapop- 
ulation are irrelevant on the global scale. A similar distinction between global and local epi- 
demics is also important when analyzing epidemics in static networks of individuals with com- 
munity structure; epidemics confined in a single community are not usually regarded as global 
epidemics [55] . 
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Figure 1: Results for the SIS model on the scale-free metapopulation network with = 200 
and (k) =6. (a) Relationship between (5^ and D for a generated network, (b) Fraction of 
infected individuals after transient, (c) Inverse participation ratio of the distribution of the 
infected individuals after transient. 
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Figure 2: Results for the SIR model on the scale-free metapopulation network with = 200 
and (k) = Q. (a) Final fraction of recovered individuals, (b) Inverse participation ratio of the 
distribution of the recovered individuals. 
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Figure 3: Results for the airport network having = 332. (a) Relationship between /3c and 
D. (b) Fraction of infected individuals in the SIS model after transient, (c) Final fraction of 
recovered individuals in the SIR model. Because of the relatively large A^, we carry out 1000 
realizations for a given combination of D and (3 to generate (a) and 500 realizations to generate 
(b) and (c). 
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Figure 4: Results for the regular random graph with = 200 and ki = (k) = 6. (a) Relation- 
ship between /3c and D for a generated network, (b) Fraction of infected individuals in the SIS 
model after transient, (c) Final fraction of recovered individuals in the SIR model. We carry 
out 5000 realizations for a given combination of D and /3 to generate (a) and 2000 realizations, 
i.e., 5 realizations for each of 400 networks, to generate (b) and (c). 
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